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NOMENCLATURE 
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Pmin 

Pmax 
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Re x 
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amplitude of the wavy surface 
Fourier component of amplitude 
Fourier-Tchebycheff expansion coefficient 
coefficient of pressure 
domain of analysis 

right-hand side of Poisson equation (See Eq. 20) 
Jacobian of transformation 
X-decomposition of function f (See Eq. 23) 

Fourier transformed variable 
transformation coefficient 
mode number 
wave number 

length scale used in squeezing transformation 

mode number 

mode cut-off number 

Tchebycheff mode number convective terms (See Eq. 11) 

Tchebycheff cut-off mode number 

pressure 

minimum pressure 

maximum pressure 

Fourier transform of pressure 

free stream pressure 

region of analysis 

Reynolds number based on length X 

time 



u 

V 

X 

X 

y 

V 
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n 

<f> 

K 

x 

V 

a 

x 
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velocity vector 

expansion function 

distance in the flow direction 

distance in transformed plane 

distance in the normal direction 

distance in the transformed plane 

boundary layer thickness 

coordinate in conformal mapped system 

phase angle 

eigenvalue 

Eigenvalue 

viscosity 

parameter in viscous stabilization 
stress tensor 

coordinate in conformally mapped system 


Superscripts* 

'v explicit convective updated variable (See Eq. 10) 

^ implicit values afterconvective update 

* pressure updated value 

- viscous update (explicit) See Eq. (25) 

— viscous stabilized values See Eq. (26) 



INTRODUCTION 


There is considerable experimental evidence (1, 2) that 
turbulent flow over wavy surfaces reduces skin friction drag size- 
ably. However, a favorable net drag reduction over wavy surfaces 
have not yet been reported mainly because of a pressure 
drag contribution in all of these experiments. In order to design 
wavy surfaces which have large skin friction drag reductions and 
insignificant pressure drag contribution, wave geometries of con- 
siderable complexities may need to be considered. The present 
investigation under NASA Contract NAS 1-15754, was motivated large- 
ly from these considerations and as an attempt to understand the 
underlying mechanism governing the details of flow changes over 
wavy surfaces through numerical simulation. 

Work on flows in wavy geometries began with the classical 
Kel vin-Helmholtz analysis of linearised inviscid flows 3 . For a 
wall displacement of the form y = a cos(kx + $) [where the assump- 
tion ka « 1 has to be made for the linearisation to hold] the 
Kelvin-Helmhotz theory predicts a pressure distribution at the wall 
of the form 

P = Pc + pU* ka cos ( kx + tt + #)• • (X) 

thus the pressure coefficient ! Cp = p - p w is 2ka at the 

— nr 1 

wal 1 . JspUoo 

The inviscid predictions overpredict the magnitude of the 
pressure distributions and have incorrect phase relations between 
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the pressure and the wall, for real flows. In order to improve the 

prediction techniques Miles^ extended the inviscid theory to account 

for shear in the mean flow. Miles' theory improves the prediction 

of pressure coefficient; however, due to the inviscid nature of the 

analysis, the pressure distribution is still 180° out of phase with 

5 

the wall. Benjamin analysed high Reynolds number laminar flows 
over wavy walls including the effect of shear through asymptotic 
theory (Ka«l). The pressure distribution at the wall, as predicted 
by his theory, has the correct phase and magnitude (for small ampli- 
tude waves in fairly high Reynolds number laminar flows). 

There have been a large number of investigations of turbulent 
flow over gentle waves. A majority of these investigations consider 
that the wall perturbation (i.e. waviness) does not modify the turbu- 
lent process as such. Davis^ has postulated that perturbation turbu- 
lent stresses be also included in analysing the problem. Zilker et 
al^ have incorporated effects of self-generated pressure gradients in 

O 

the length-scale models. Recently Markatos has attempted the solution 
of wavy wall turbulent flow through a two equation turbulence model by 
finite difference techniques. While Markatos 1 numerical results are 
at best first order accurate in space, they go a long way towards 
providing understanding of the nature of flow over wavy walls. 

The present work extends previous investigations in several 
ways. The full time-dependant two dimensional Navi er Stokes 
equations are solved using spectral methods (to achieve high spatial 
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accuracy) and higher order time splitting methods together with con- 
formal mapping methods to allow simulation of flow over steep 
waves (ka ^ 0 (1)) Two basic computer codes have been developed 
under the present contract. The two specific codes described in 
this report are: 

1) Code 1 - Periodic inflow-outflow conditions are used 

in the lateral direction for walls with spatial regularity. 

A spectral Tchebycheff-Fourier mode expansion is used 
(spectral Tchebycheff--in the normal direction of the flow, 
and Fourier— spectral in the direction of periodicity) to 
obtain infinite order accuracy in space. A second order 
accurate time scheme (Adams-Bashforth with Richardson 
extrapolation) is adopted in this code to obtain accurate 
time transients. 

2) Code 2 - This computer code uses prescribed inflow-outflow 
conditions; a two dimensional Tchebycheff expansion scheme 

is used to guarantee high spectral accuracy. Again, second 
order time accuracies are obtained using Richardson extrap- 
olation schemes. 

In Section II the methods of solution adopted for the two codes 
are described. Results for laminar flows at various Reynolds numbers 
are presented in Section III. In Section IV, results for nonsymnetri c 
waves are described. Finally, in Section V application of these codes 
to other flow problems and directions of future work are delineated. 
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II. METHODS OF SOLUTION 
The two dimensional Navier Stokes equations are 

3u_ + u. Vu = - Vp + V.t m 

3t ' 

V.U = 0 (2) 

where the stress tensor x may include both laminar stresses 
and turbulent stresses. 

II. A METHOD OF SOLUTION FOR CODE I 

For Code I the domain of analysis and the boundary condi- 
tions appropriate are 

- D: 0 < X < 2tt , f(x,t) <.y 1“ (3) 

above the wall, y = f(x,t) 

and u(x,y,t) = u(x+2ir J y, t) 

p(x,y,t) = p(x+2iT,y,t) 

"r (x,y,t) = r (x+2 7r,y,t) 

The test results reported (using Code I) i.n this report are for the 

case where the wall is stationary, i.e. 

% 

f (x >t) = Z =1 a n i cos(k n x + <|> n ) ( 5 ) 

where m is the cut off mode number in the summation. A conformal 

Q 

mapping technique is used to transform the domain (3) into the 
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region* 


( 6 ) 


0 £ <; £ 2 tt , 0£n£°° 

When the slope of the wall |ka| £ 1 the conformal mapping 
coefficients are accurately generated in only N log N operations 
where N is the resolution along x(t). Thus time-dependant 
geometries can be very efficiently handled through this mapping 
technique. A further stretching (non-conformal ) transformation of 
h is carried out in order to implement the spectral method so 
that the domain (6) maps into the finite region. 


R: 0 < X < 2u , -1 < Y < + 1 


(7) 


by the transformation 


x = d 


, Y = 2n - 1 
n+L 


or 


X = c 


-n/L 

, Y = 1-2 e 


( 8 ) 


where L is some suitable mapping scale. 

A plot of the mapped grid used in the computation is 
shown in Fig. 1. The wall is a sine wave (m - 1 ,i <!> n = * ). The 
wall slope ka = 0.1963. Only a portion of the computational grid 
is shown in Fig. 1 (so as to keep x, y scales in plotting nearly 
equal). The vertical lines in Fig. 1 are lines of constant C , 
and the horizontal lines are lines of constant ( n : Y ) The 
orthonormality of the grid lines is quite apparent from the figure. 
Equations (1) — (2) are solved numerically by a fractional step 
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procedure. The procedure adopted is the spectral collocation tech- 
nique. In this method, all dependant variables are expanded in 
Fourier Tchebycheff series as 

c/v » N imX 

F(X,Y,t) = Z Z A(n,m,t)e T (Y) (9) 

n=0 |m| < m c n K ' 

where ^" n M = cos(n cos* Y) is the Tchebycheff 

polynomial of degree n. The fractional steps used in this pro- 
cedure are: . 

(i ) Convective update 
(a) Explicit step 

An intermediate field is found such that 

J = U n (x ' y>t) *3 it N(x,y,t„) (10) 

where 


N(x,y,t) - - (u - U m (Y)).Vu (n 

Urn in Eq (11) is a mean profile (of the dominant velocity component 
in X direction) 


(b) Implicit convective step 

y v 

u + At Urn * Vu = u 


(12) 


The application of the implicit step releases time step restrictions 

of Courant-Friedrichs-Lavy type (At < a* ) 

U 
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( i i ) Pressure update 


% 

Oi 


u = u - vp 


(13) 


A 

U 


= 0 


(14) 


with the boundary conditions 

/v 

n • u =0 (at n - 0) 


(15) 


from (13), 


V p = v • u 


(16) 


(for a stationary wall) 


(17) 


with 3£ - 0 

an 

The boundary conditions in c directions are already satisfied by the 


assumption of periodicity of the flow. 

( i i b ) Fast poisson solver 

Solution of (16) with (17) is implemented through a fast 
Poisson algorithym. The conformality of mapping between (x, y) to 
(? , n) provides 


v 2 P = |F'l* (ai+ifP) 

3X* arf 7 

The right hand side of (18) can be written as 

v2p = (F') 2 (3jp + g(Y) 8 (g(Y)3£)) 

? X 3V 3Y 


where g(Y) = 3Y 
9n 


(18) 


( 19 ) 
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Equation (16), therefore, modifies to 


ifP + g(V)3 (g(Y)9j3) = f(X,Y) (20) 

aF 3V 3Y 

In order to solve Eq (20) efficiently, we introduce the eigen- 
value problem. 


g(Y) 9 (g(Y) 9uJ = Xu (21) 

9Y 9Y 

subject to the homogeneous boundary conditions, 


^-0 9 Y - -1 ( 22 ) 

3Y 


The eigen values and the eigenfunctions associated with this problem 
(21) are found in a pre-processing step. The function u(X, Y) can 
therefore be expanded as. 


u(X,Y) = Z V k (X) u k (Y) 
k 


(23) 


and similarly 

f (X ,Y) = T. F k (X) U k (Y) 

Thus Eq (20) can be written as 

,2, I ((1 -l) 2 + X k ) p = F 

|1|= 0 k - 0 (24 ) 


where the til da stands for fourier transformed variables. Solution 
of the poisson equation is thus accomplished through further matrix 
operations of Eq (24). 
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( i i i ) Viscous update 


u = u + At V • T n ( 25 ) 

with the boundary condition 
lT = 0 @ n = 0 


This Euler step is only first order accurate but will be improved 
upon later. The time step restrictions for the numerical stability 
of (25) are particularly severe, especially in the case of turbu- 
lent flows so the next fractional step is designed to stabilize the 
scheme. 

(iv) Viscous stabilization 



wi th 

IT = 0 


At' ) V 2 U- = u - ( oAt v ) v 2 u n 

|F’ I 2 


@ T) = 0 


(26) 


Here ' a is an optimising parameter for stability and accuracy. 

The stabilization step is implemented along the lines of the poisson 
algorithm (i.e. via eigen function expansions) 

(v) Local extrapolation 

For time dependant problems requiring accurate time tran- 
sients, a local extrapolation procedure (10) is used to achieve 
global second-order accuracy in time. 
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II. B METHOD OF SOLUTION FOR CODE II 


The basic algorithm for Code II contains the five fractional 
steps discussed above. There are qualitative differences in Code I 
and II; these differences arise because the boundary conditions at 
X - 0, 2tt are different. Because Tchebycheff polynomials are 
used in X direction also, the equation (24) for poisson equation for 
Code II is of the form 


Nc N 2 c 

2 Z + X n ) Pmn = Fmn 

m=0 n=0 


(27) 


where < are the eigen values in X direction and X are the 
m n 

eigen values in Y direction; j , F^- are the decomposed modal 

components of p(X,Y) and f(X,Y) in the eigen space (m, n). 


In the next section results are reported from these codes 

using, 

(i) up to 32 Fourier modes and 64 Tchebycheff modes for’ Code I. 

(ii) 16 Tchebycheff and 32 Tchebycheff modes for Code II. Typically 
the computer codes run at 2 ms per retained mode per time, step on 
CDC 6600 computer. 
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III. LAMINAR TEST RESULTS 


A large number of laminar flows over wavy walls were done 
using Codes I & II. Results from three representative runs will 
be described here. 

Run 1 - The flow situation here is over a very gentle wave 

(wave surface y = a cos(kx) where a = 2.54 u , k = 1.2566 mm”' , 

-3 5 

ka = 3.19 x 10 ) at a Reynolds number Re = 4 x 10 with a starting 

length at entrance X = 0.287m, so the boundary layer thickness is 0.3 cm. 

This experiment is appropriate for comparison with Kelvin-Helmholz 

theory and Benjamin's theory (which requires ka«l. Re »1). For this 

A 

_3 

flow, the Kelvin-Helmholz inviscid analysis predicts C = 2ka = 6.3 x 10 . 

r 

The pressure distribution predicted for this flow by this theory is 180° 

out of phase with the wall. On the other hand, Benjamin's theory pre- 

-4 

diets a phase shift of 217° and C p = 1.8 x 10 . The results obtained 

using Code I indicates excellent agreement with Benjamin's theory. The 

phase of the pressure distribution was found to be 214° ± 5° behind the 

surface displacement and the pressure coefficient was found to be 
-4 

Cp = 1.742 x 10 at the wall. In Fig. 2, the contours of pressure at 
steady state are p 1 otte d . The rate of decay of the pressure distribu- 
tion away from the wall in this calculation is not in agreement with the 
predictions of Benjamin's theory. Close to the wall the pressure decays 
very slowly ;at ky = .78 tne pressure has decayed only by a factor 

0.8238, so p 'v — y — , whereas the theory predicts decay of the pressure 
p e”^. At intermediate distances away from the wall, the computer 
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calculations show that the pressure decay returns to e~^ gradually. 
In Fig. 3, contours of velocity in the y direction are plotted. It 
is apparent from the figures that the y - velocity distribution peaks 
at a small distance from the wall and has phase lag of 63® with the 
surface. 

Run 2. - The second run is similar to the first run except 
that the amplitude of the surface displacement is increased to 
ka = 5.5 x 10~ 2 . Kelvin-Helmholz' theory predicts Cp = 0.11 for this 
case indicating that the wave amplitude is much too large to be ana- 
lysable through pure inviscid theory. The Benjamin' theory predicts 
a phase shift of 217° with Cp = 3.34 x 10~ 2 . Computer calculations 
indicate pressure coefficient Cp = 3.3 x 10 and a phase lag of 
215 - 6°. Again the rate of pressure decay near the wall disagrees 
with theoretical models. 

Run 3 . - The third run involves similation of experiments by 
Kachanov et. al.^ Here k = 25mm~^ while ka = (5.3 ± .01) x 10 -2 . 

The Reynolds number for the flow is Rex = 1.1 x 10 5 . U = 5.3 m/sec. 

CO 

The boundary layer thickness at the start of the computational 
domain is .6 cm where X = 3.06 cm. In Fig. 4, the contours of the 
x-component of the velocity are plotted at a steady state. From 
Fig. 4, it is seen that there is a thickening and thinning of the 
boundary layer in the trough and crest regions of the wave. In 
Fig. 5, contours of the y-velocity are plotted; as before, the 
maxima and minima lag in phase by about 63°. In Fig. 6, contours 
of the pressure in the unmapped x-y domain are plotted from the 


13 



figure, the maximum pressure is found to be considerably smaller 

than the minimum pressure so the disturbance due to the wall is no 

-3 -3 

longer linearizable; in fact, p = 5.2 x 10 while p m ^„ = 4.5 x 10 . 

Comparison of the fluctuation amplitude Au at the height of unper- 
turbed boundary layer thickness between the experiment and our 
calculations indicates excellent agreement. 
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IV. NON SYMMETRIC WAVES 


The flow over nonsymmetrical periodic waves can also be 
computed using the present code. Here we report on calculations 
for two such nonsymmetrical waves. The wall amplitude and wave- 
length are the same as that in the Kachanov et al^ experiment. 

Also, the flow conditions are kept the same as in the latter experi- 
ment in order to see whether the flow over nonsymmetrical waves 
can be beneficial for drag reduction. In Fig. 7 the grid lines 
are shown for one such nonsymmetrical wave. The conformal nature 
of the transformation is not apparent in this grid plot because 
the x, y length scales are different. In Fig. 8, contours of the 
x - velocity are plotted. As before, we see a thickening and 
thinning of the boundary layer thickness corresponding to decreasing 
and increasing wall skin friction. The noisy contour near the 
high is due to the sensitivity of the contour to small discretiza- 
tion errors. In Fig. 9, contours of the y-velocity are plotted. 
There is a marked change in the phase shifts of the peaks and 
valleys of the velocity relative to the wave surface. In Fig. 10, 
pressure contours are plotted. The magnitude of C is about the 

r 

same as in Kachanov et al experiment. However, the pressure 
distribution has a maximum almost 180° ahead of the wave peak and 
a minimum about 80° ahead of the trough. Obviously, the skin fric- 
tion and pressure distributions for this flow are quite different 
from those of the symmetrical wave surface. Nevertheless, there 
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are only small changes in both the integrated friction drag and 
pressure drag from the Kachanov et al experiment, which we estimate 
has a net drag reduction of about 5 % compared to a flat plate. 

In Fig. 11, the grid system for the second nonsynmetri c 
wave is plotted. Again, the flow conditions as well as the wall 
amplitude and wavelength are the same as in the Kachanov et al 
experiment. In Fig. 12, contours of the x-velocity distribution 
at a steady state are plotted. The thickening and thinning of the 
boundary layer profiles is again apparent. In Fig. 13, the y-velocity 
is contoured. It is apparent that the highs and lows of the y-velocity 
are markedly different in both magnitude and location. Because 
of the long duration of favorable pressure gradient, the activity 
is quite low (small negative y-velocity) in the cusp regions. The 
steep rise of adverse pressure gradient causes increased activity on 
the lee-side of the wave. 

In Fig. 14, the pressure distribution for this wave is 
presented. Again, the amplitude and phase of the pressure distri- 
bution is markedly changed. However, the integrated frictional 
drag and pressure drag show no significant change from the Kachanov 
et al , experiment. 

Flow over Bumps 

In Fig. 15, 16, 17 the flow calculations over a finite 
bump are reported. The bump is Gaussian and Fig. 15 shows the 
computational grid system. In Fig. 16, the pressure contours are 
plotted. Fig. 17 shows the x-velocity distribution over the bump 
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V. FUTURE DIRECTION FOR WORK 
APPLICATIONS OF THE NAVIER-STOKES CODES TO OTHER PROBLEMS 

The computer algorithms developed in this report have 

been used to study turbulent flow over wavy walls. The experi- 

1 ? 

ments chosen to study have been that of Kendall and Si gal . 
Currently, the effects of turbulence are solved outside the code 
itself via a boundary layer program which uses multi -equation 
turbulence modelling. The mechanics of such analysis are described 
below. 

(i) The boundary layer program requires the self generated 
pressure distribution which is obtained from the Navier Stokes Code. 

(ii) The Navier Stokes Code uses the lumped eddy stresses 
from the boundary layer program. Fig. 18 shows the pressure distri- 
bution obtained by the Navier Stokes Code using this technique over 

a Kendall surface. The experimental value of pressure coefficient 
for this problem is Cp = 0.1; our calculations show a pressure 
coefficient of Cp = 0.13. The phase of the pressure distribution 
obtained by our analysis compares extremely well with the experi- 
ment of Kendall. The disagreement in the amplitude of the pressure 
distribution indicates that, turbulence modelling be done directly 
within the Navier Stokes Code itself rather than through the 
boundary layer program. Basically a time transient code such as 
ours changes the mean flow profiles with time and hence the 
turbulent stresses have to change accordingly. The interaction 
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technique described above is inadequate to handle this situation. 

One of the chief dilemmas concerning incorporation of 
turbulence modelling in a Navier Stokes Code such as this revolves 
around the inadequate experimental data over which such models may 
be compared with. 

The inflow-outflow code (Code II) will be an excellent 
start off point for such investigations. Unfortunately, contractural 
constraints of finite time objectives often leave such work in 
jeopardy. 

The inflow- outflow code can also be used to analyse con- 
verted eddies within the computational domains. Presently such 
investigations are being conducted at LaRC. Due to the extreme 
sensitivity of the code to smooth initial data, such experiments 
can not be reported at the present time. 


This work was supported by NASA Langley Research Center 


under contract NASI -15754. 
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FIGURE TITLES 

Figure 1. Plot of Mapped grid system used In computation. 

The wave surface is a sine wave (Ka = 0.1963). 

The orthonomality of grid lines is apparent in 
this plot. 

Figure 2. Contours of pressure at a steady state for Run I 
• -3 

(ka = 3.19 x 10, flow is laminar at Re^ » 4 x 10°) 

Figure 3. Contour of y-velocity distribution for Run I. 

Figure 4. Contours of x-velocity at a steady state for the 
experiment of Kachanov et al". (Run III). 

Figure 5. Contour of Y-velocity for Run III. 

i 

Figure 6. Contours of pressure for Run III. 

Figure 7. Computational grid system for nonsymmetri c wave 
(Laminar flow Re x - 1.1 x 10 5 ). 

Figure 8. Contours of x-velocity over nonsymmetric 
wave. 

Figure 9. Contours of y-velocity over a nonsymmetric wave. 
Figure 10. Contours of pressure for flow over nonsymmetric 


wave. 



Figure 11 


Figure 12 


Figure 13.. 


Figure 14. 


Figure 15. 


Figure 16. 


Figure 17. 


Figure 18. 


. Grid system for the second nonsymmetri c wave 
(Laminar flow, Re x = 1.1 x 10 5 ) . 

• Contours of x -velocity at steady state for the 
second nonsymmetric wave. 

' Contours of y-velocity at steady state for the 
second nonsymmetric wave. 

Contours of pressure at a steady state for the 
second nonsymmetric wave. 

Computational grid system for flow over a periodic 
ga us si an bump. (Laminar flow Re y =1.1 x 10 5 ) V 

Contours of pressure at steady state for flow V 
over a periodic bump. ^ ‘ ; \ 

Contours of x-velocity for flow over a periodic 
bump. 

Contours of pressure at steady state for turbulent 
flow over wavy surface corresponding to the experi- 
ment of Kendall 
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The wave surface Is a sine wave (Ka = 0.1963). 
The orthonomall ty of grid lines Is apparent In 
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Figure 2. Contours of pressure at a steady state for Run I 

, -3 c 

(ka - 3.19 x 10, flow Is laminar at Re v = 4 x 10°). 


V - VELOCITY PROFILE 



Figure 3. Contour of y-veloclty 


for Run I 
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Figure 6. Contours of pressure for Run III. 
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Figure 7. Computational grid system for nonsymmetrl c wave 
(Laminar flow Re x = 1.1 x 10 5 ). 
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Contours of pressure for flow over nonsymmetrlc 
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Figure 11. Grid system for the second non symmetric wave 2 * x , > 


(Laminar flow, Re x = 1.1 x 10 5 ). 
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Figure 13. Contours of y-veloclty at steady state for the 
second nonsymmetrl c wave. 
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Figure 15. Computational grid system for flow over a periodic 
gausslan bump. (Laminar flow Re v = 1.1 x 10 5 ) 



Figure 16. Contours of pressure at steady 
over a periodic bump. 
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Figure 17. Contours of x-veloclty for flow over a periodic j 
bump. ! 
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Figure 18 . Contours of pressure at steady state for turbulent 
flow over wavy surface corresponding to the experl 
ment of Kendall! 1 ) 
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